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Abstract 

The distribution of net electric charge in graphene is investigated, using both a constitutive 
atomic charge-dipole interaction model and an approximate analytical solution to Laplace's equa- 
tion. We demonstrate a strong size dependence of the charge distributions in finite-size, infinitely- 
long and multi-layered rectangular graphene sheets, respectively. We found that the charge den- 
sity can be naturally enhanced up to 13 times at graphene's geometry edges. This edge charge 
enhancement effect becomes more significant when the length or the width of graphene increases. 
The charge enhancement ratio is found to follow a linear relationship with the number of layers. 
These results can be used to understand the newly experimentally observed electron emission, 
charge impurity and chemical doping phenomena in 2-dimension nanostructure. 
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I. INTRODUCTION 



Graphene is an hexagonally coordinated semi-metallic or zero-gap-semiconducting carbon 
layer that crystallizes in sp 2 structure. Its unusual band structure [H |2] leads to unique 
electronic properties such as high charge carrier mobility [31 HI [5] and tunable electron 
band gap by external electromagnetic fields j6j [7J, El El QUI E] or by chemical doping. [121 02] 
These make graphene an ideal candidate for a number of applications in nanoelectronics. [H] 
Recently, it has been shown by electrostatic force microscopy experiments that electric 
charges can be stored in carbon nanotubes for a period of time. [151 HS1 Ej It would not be 
surprising if graphene exhibits similar properties, because of the large similarity in its atomic 
structures to that of nanotubes. [T8] The investigation on electronic properties of graphene 
doped with net charges is of importance for understanding the charge impurities [191 EH EI] 
or chemical doping [121 03] in graphene, and may as well lead to further applications in 
charge storages, [22] sensors [23J and field-emission devices [2J1 ES] based on graphene. 

The question about the exact solution to the charge distribution in an arbitrary-shaped 
conducting surface has been raised since the theoretical basis of electromagnetism was estab- 
lished by J. C. Maxwell in 1873. [26] For graphene, the electrostatic problem of finding the 
charge distribution on the surface of a conducting, rectangular plate requires the solution to 
Laplace's equation in the three-dimensional geometry. Unfortunately, even such a seemingly 
simple problem has no exact solution: The number of exact solutions in three-dimensions is 
still limited to specialized geometries that exhibit a great deal of symmetry [271 EE1 El] • In this 
paper, the charge distribution was first calculated using a numerical atomic charge-dipole 
model. Comparisons are then made between the result of a simple, one-term variational 
approximation and a numerical solution, which is found to give a more accurate description 
of the edge atoms. 

The remainder of this paper is organized as follows. In Sec. (|n]) we describe first our 
numerical method and then we demonstrate a variational approach to solve Laplace's equa- 
tion and calculate the charge distribution on a conducting rectangular plate. The results of 



both numerical and analytical models are compared and discussed in Sec. (III). We draw 



a conclusion in Sec. (IV). Mathematical details of the one-term variational approximation 
are summarized in the Appendix. 
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II. THEORY 



A. Numerical method 

In this work, the charge density was first calculated using a computational atomic charge- 
dipole model, [3n] in which each atom is modeled as a polarizable sphere with a net charge 
q and an induced dipole p. The equilibrium distribution of the charges and the induced 
dipoles in graphene is determined by minimizing the total electrostatic potential energy of 
the system 

N N N N 

jjelec = j2 qi ( Xi + Vi) - E Pi ■ Ei + \ E E QiT^ q q, 

i=l i=l i=l j=l /i\ 

N N N N I 1 / 

i=lj=l i=l j=l 

where x is the electron affinity, and where V and E are the external potential and electric 
field, respectively. The electrostatic interacting tensors are T and T, which are regularized 
by a Gaussian distribution [31] . Requiring the derivatives of JJ elec with respect to q and p 
to be zero, Eq. (jl]) yields 

N N 

Y.T;%®p 3 + T.T l p l q q 3 = -E i 

3=1 3=1 
N N 

E ni q ■ Pj + E r q l q q 3 + a = —(xi + Vi) 

i=i 3=1 (2) 

N 

E Qi = Q tot 

Vi = 1,...,N 

where A is the Lagrange multiplier which is related to the chemical potential of the molecule. 

This numerical model is an extension of the atomic dipole interaction theories [321 [331 El] ■ 
It has recently been parameterized for sp 2 carbon systems [3S1 EE] and validated through 
electrostatic force microscopy charge injection experiments. [3T] The graphene structure is 
relaxed by means of an atomic energy-optimization method using the adaptive interatomic 
reactive empirical bond order (AIREBO) potential function. [38] which has specially been 
parameterized for sp 2 hydrocarbon systems. [391 EH HH H21 H3] 
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B. Analytic approximation 



A simple analytical approximation for the solution to Laplace's equation is now con- 
structed. In three physical dimensions, consider an infinitesimally thin conducting surface 
of rectangular shape L x -by-L y that is held to the potential Vq, with the reference or zero 
potential at infinity. The rectangle lies in the z = plane and is the region \x\ < L x /2, 
\y\ < Ly/2. We seek a solution to Laplace's equation 

d 2 d 2 d 2 \ 

+ 7^ + 7^ hKw)=0 (3) 



(4) 



dx 2 dy 2 dz 2 / 
subject to the boundary conditions 

il>{x,y,0) = V (\x\ < L x /2, \y\ < L y /2) 
ip(x, y, z) — > as a/x 2 + y 2 + z 2 — > oo 

Based on experience with the exact solution for the two-dimensional strip problem, 
we propose a one-term approximation for the (scaled) charge density 

{—k 
— . — , \x\ < L x /2 and \y\ < L y /2 

y/1 - (2x/L x ) VI - WJLJ 2 (5) 
0, \x\ > L x /2 or \y\ > L y /2 

where the single unknown is the constant k. This functional form is simply a factorization 
of the rectangular distribution into two edge-condition forms that apply to the classical two- 
dimensional strip problem. Its suitability is not a priori known, but ultimate comparison 
with independent moment method [HZ] calculations shows that it is quite respectable. The 
one-term Galerkin procedure in the Appendix produces an approximation for the constant 

2V 



(6) 



in terms of an integral function I(L X /L y ) of Eq. (25) which is evaluated numerically. The 
resultant surface charge density is then 

<*(x,y) d . 2k 

— — 2— ib(x, y, CT) = . — . (7) 

^ V 7 ! " (2x/L x ) VI " ^y/L y ) 2 

Integration of Eq. ([7]) gives the total charge 

Lx/I Ly/2 



induced on the rectangular plate when held to the potential Vq. 

Here Fig. [T] shows the charge enhancement effect in a rectangular plate exhibited by 
Eq. ([7]). In addition to the surface charge distribution, it is instructive to also consider the 
capacitance C = Q/Vq, a global electrostatic characterization of the plate. In the present 
notation, the capacitance per unit width is 

- = — (9) 

In the limiting case of a square, L x = L y and the numerics give 1(1) = 2.176 and the 
subsequent normalized capacitance is (in SI units) C/L y = 40.159 (pF/m), about 1.6% 
lower than the value 40.811 (pF/m) that has been published by Goto et al. |44j 
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III. RESULTS AND DISCUSSIONS 



A. Finite-size graphene 

We start with mono-layered graphene which contains a quantity of net electric charge. 
Fig. |2] (a) shows the charge distribution in a finite-size graphene sheet, calculated using the 



numerical method (Sec. (II A)). As expected, typical charge and electric field enhancement 
effects [37] were observed near the corners and the edges. To compare these results with 
those predicted by Eq. (|7]), we show the atomic charge density in a quarter of the graphene 
sheet in Fig. [2](b). Comparison between the contours of charge density in Fig. [ljand Fig. [2] 
(b) shows a general agreement between the analytical model and the numerical calculation. 
The lack of smooth trend of the contours in Fig. [2] (b) is due to the different configurations 
of atoms at armchair and zigzag edges and the discontinuity of our discrete model. It is 
found that the Eq. ([7]) becomes no more accurate for the atoms near the edges, because the 
assumption of infinitely sharp edges in the analytical solution is not valid in atomic scale. 
After our simulation results, the charge density (at the corner) can be enhanced up to 14 
times over that in the center of graphene. We note that the charge density can also slightly 
vary if the graphene is substitutionally doped by net charges (e.g. some electrons per 100 
atoms), since the graphene band structure can be significantly modified in highly charged 
graphene. jlH] Moreover, these results also imply that the influence of induced dipoles on the 
charge distribution is minor. 

We next consider the size effect on the charge distribution. Our results show a strong 
dependence of the charge enhancement on the graphene length L x (Fig. [3J. We can see that 
the fluctuation of the charge profile due to the edge state becomes less significant with the 
increasing length, and that the charge enhancement at the edges is stronger for the longer 
sheets. This behavior is similar to that in carbon nanotubes. (13] 

We now define the maximal charge enhancement ratio 7 maa; as the maximal atomic charge 
density normalized to the average over the whole graphene. It is one of the key parameters 
for the potential applications of graphene in field emission devices |46j . To show a complete 
correlation between jmax and the dimensions of graphene, we plot numerical calculation data 
of 'Jmax m Fig. [4] as a function of L x and L y . We can see that 7 ma x increases with either 
L x or L y , and that 7 ma x gets higher for L x « L y . It is found that the increase of charge 
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enhance ratio is more rapid in small size graphene sheets. 



B. Infinitely-long graphene 

To generalize our results for the graphene size usually used in experiments, it is partic- 
ularly interesting to investigate the charge distribution in infinitely long graphene sheets 
or strips. The question of charge distribution becomes a 2D standard problem in elliptical 
coordinates. [28J The exact solution to Laplace's equation for the strip conductor of width 
L x held to the potential Vo gives the surface charge density 

(Tlx) = . , (10) 

L x In (Ap /L x ) y/1 - (2x/L x ) 2 

where po is the (large) radial distance where the potential vanishes. This is the "ground" or 
reference surface that cannot be at infinity because of the logarithmic singularity in the 2D 
potential problem. Integration of the surface charge density yields the capacitance per unit 
length 

Ci = ln(2 2 po/l) ' 

Here we show the normalized charge profile in a graphene sheet with applied periodic con- 
dition in width (y) direction in Fig. [5j which highlights a good agreement obtained between 
the numerical calculations and the analytical solution. Also, compared to Fig. |3j it can 
be seen that the fluctuation on the charge density curves due to the open-edge effects [H] 
completely disappeared in the infinite-long graphene. [I8l H9] 

To demonstrate the size effect, we now consider the maximal charge enhancement ratio 
Imax as a function of L x . It is shown in Fig. [6] that ^ m ax increases with L x , the shape of the 
Imax — vs — L x curve is as predicted by the logarithm relationship proposed by Keblinski et 
al. jl5]. Since the surface electric potential is a constant, this length dependence directly 
leads to the rise of the capacitance C of the edge atoms. Also, we plot the normalized C in 



Fig. |6j It can be seen that C also increases with L x , as predicted by Eq. (11). However, 
it is found that the capacitance per unit length decreases with L x . This behavior indicates 
an electric potential loss when the graphene sheet gets shorter, i.e., the longer graphene 
exhibits a higher surface potential to keep the net charges inside. 
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C. Mult i- layered graphene 



Multi-layered graphene sheets have also been studied in this work. They have been 
reported to have interesting electronic properties. [SUl EE] We investigated multi-layered 
graphene in which each layer contains the same quantity of net charges. Due to the electro- 
static interaction between the layers, the charge profile in multi-layered graphene cannot be 
assumed to be a simple addition of that of a single layer. In Fig. [7], we plot the profile of 
net charges in a 4-layered graphene to compare with that in a single-layered one with the 
same size. It can be seen that the charge enhancement is more significant in multi-layered 
graphene. The mechanism of this behavior relies on the electrostatic repulsive effects be- 
tween the carbon layers, which can effectively reduce the capacitance of the atoms at the 
center of graphene and can however enhance that of the edge atoms. [52] 

We now estimate the charge enhancement ratio in multi-layered graphene as a function 
of the number of layer N (Fig. [8]). Our results showed that the charge enhancement 
effect becomes more significant when the number of layer increases, due to the electrostatic 
repulsive effect between the layers. We found that 7 max roughly follows a linear relationship 
with N. Note that the interlayer distance can be slightly changed in case of a very high 
charge density. [53] while we use a constant distance 0.34 nm in this work. [54] 



8 



IV. CONCLUSIONS 



In conclusion, we have studied the charge distribution in graphene, using both a numerical 
method and an analytical solution to the classical problem of net charges in a rectangular 
conducting plate in three physical dimensions. It is found that the analytical model gives 
a good description to the charge distribution, except for the very edge atoms because of 
the atomic edge state. We demonstrated how the charge distribution changes with the 
increasing graphene length, width and number of layer. The length dependence of the 
charge enhancement effect was found to be related to the rise of graphene capacitance with 
increasing size. The charge enhancement ratio is found to follow a linear relationship with the 
number of layers in multi-layered graphene. These results developed a better understanding 
of the state of net charge in graphene, which can be useful for the potential applications of 
graphene in charge storages, sensors or field-emission devices. 
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Appendix: Variational approximation for the charged rectangular plate 



The desired slotion to Laplace's equation Eq. (|3| subject to the symmetric boundary 
conditions Eq. Q is even in x and in y (and in z), so Fourier cosine integrals 



OO OQ 



ip(x,y,z) = J da J df3F{a, (3) cos(ax) cos(f3y)e 7|z| (12) 
o o 

are appropriate with a 2 + (3 2 — 7 2 = 0. The non-analytic modulus \z\ accommodates the 
discontinuity in normal (z) derivative 

OO OO 

— ip(x, y, 0*) = =F J da J dj3 ^F^a, (3) cos(ax) cos(j3y) (13) 
o o 

across the surface charge. For a and (3 real, choose 7 = + \Ja 2 + (3 2 . Using Eq. (j5j), Fourier 
inversion of Eq. ( 13 ) is 

^--^ k hv^hv^- (l4) 



7T 



The required integral is a Bessel function 

L x /2 



dx - COS(ax) = ^J (aL x /2) . (15) 



Hence the spectral function can be expressed as 



F(a 3) = kL x L y J (aL x /2)J <J3L y /2) 
4 ^a 2 + P 2 



whereupon the forced boundary condition in Eq. Q becomes 



00 00 



kL x L y f f J (aL x /2)J (f3L y /2) \ \x\ < L x /2 

1 da I dp cos(ax) cos(py) = V \ • (17) 

V® 2 + P 2 \\y\< L y /2 





Application of the Galerkin inner-product operator 

L x j 2 Ly 1 2 

dx f dy 







y/1 - {2x/L x f J - (2y/L y y 



to each side produces 

00 00 

kLr r . r , j2(u)j2(v 



(18) 



du I dv -wo^ = = v (19) 
' V^ 2 + (W 
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where u = aL x /2 and v = (3L y /2. With /i = L x /L y , we require the double integral 

oo oo 

I(H) = j du j dv Jl{u)Jl{v) [u 2 + (H 2 ] _1/2 • (20) 



Note that /(l//i) = fJ>I(n)- Using formulae from Watson's treatise|55j, the square of the 
Bessel function is expressed SIS db double integral 

2 tt/2 tt/2 

Jq(m) = ^— ^ y dx j d£ cos(2wcosxcos£) (21) 
o o 



of trigonometric functions, so that Eq. (20 ) temporarily balloons out into the six-fold integral 

tt/2 71-/2 tt/2 tt/2 OO oo 

r f cos(2wcosxcos£) cos( 

""' | 1/2 



m = ( -V [dx [difdy [d V [du [ d|; cos(2 M cos x cos Q cos(2r, cosy cos (22) 



7T 



[m 2 + (H 2 ] 





Formula 3.754(2) from Gradshteyn and Ryzhik[56] reveals that the it-integral is a modified 



Bessel function 



oo 

/cos(2wcosxcos£) 
du Y/2~ = -^o(2/if cos x cos 4) (23) 







[u 2 + {fiv) 2 \ 

and then their formula 6.671(14) eliminates the Bessel function in the v-integral 



oo 

dv Kq(2[iv cosxcos£) = — [(/i cos x cos£) 2 + (cosy cos 77) 2 ] . (24) 



Finally, the original double integral Eq. (20) of oscillatory functions over infinite ranges is 
converted to a four-fold integral 



of slowly-varying functions over finite ranges. A specially designed Gaussian quadrature 
routine gives numerical values with little fuss. Any one of the four integrals is a complete 
elliptic integral, but a vectorized four-dimensional quadrature is competitive in execution 
time and accuracy with a thee-dimensional quadrature that requires evaluation of the elliptic 
function. The desired constant k in Eq. (19) now takes the form in Eq. (|6]). 
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Figures 




FIG. 1: Surface charge density (a) normalized to that in the plate center. 
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FIG. 2: (Color online) (a) Representative atomic diagram of net electric charges density in a 
graphene sheet (L x x L y = 5 x 8 nm) . The color scale of atom is proportional to the charge density. 
The dark arrows stand for the electric fields induced by the net charges. The armchair edge is 
along the x axis, and the zigzag one is in y direction, (b) Atomic charge density normalized to 
that in the graphene center. 
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Normalized position along xaxis 



FIG. 3: Charge profile along the x axis in 5 graphene sheets {L y = 0.23 nm) with different 
lengths L x . The charge density is normalized to that in the center of graphene. The x position 
is normalized to L x . The symbols represent the simulation data and the dot- line curve stands for 
the analytical solution from Eq. (17]). 
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FIG. 4: Maximal charge enhancement factor ^ m ax vs. the graphene size. L x and L y stand for the 
dimensions of graphene along the armchair and the zigzag edges, respectively. 
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FIG. 5: Profile of charge density (normalized to the average) along the x axis in a charged graphene 
sheet {L x = 50 nm). The periodic condition is applied along y axis hence L y — > oo. The symbols 
represent the simulation data and the curve stands for the classical solution from Ref. [28] . 
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FIG. 6: Maximal charge enhancement factor j m ax and normalized capacitance C of graphene 
sheets of various lengths L x (L y — ► oo). 
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FIG. 7: Profile of charge density (normalized to the average) along the x axis in a multi-layered 
(filled symbols) and a single-layered (empty symbols) graphene sheets {L x = 20 nm, L y — > oo). 
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FIG. 8: Maximal charge enhancement factor 7 max vs. the number of layers N in graphene sheets 
{L x = 20 nm, L y — > oo) 
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